%% First load data arrays for analysis
clear;
fnames = dir("monod*.mat");
for k=1:length(fnames)
    load(fnames(k).name);
end


clc;
close all;

%% Establish concentrations to evaluate
ivec = [1:6]'; % Only use first 4 concentrations

%% Define molar mass
% Define approximate molar mass of components
molmass_casamino = 135.5; % g/L

% Total carbon molar mass
molmass_NZCasePlus = 125/molmass_casamino;

tempmat = monod_mat_NZCasePlus_37C;
concentration = tempmat(ivec,1)/100*molmass_NZCasePlus*1000; % Convert to mM

%% 37C NZCasePlus
tempmat = monod_mat_NZCasePlus_25C;
smat = size(tempmat);


grmax_25C = nanmean(tempmat(ivec,2:4),2);
grmax_err_25C = nanstd(tempmat(ivec,2:4),0,2);

%% 30C NZCasePlus
tempmat = monod_mat_NZCasePlus_30C;
smat = size(tempmat);


grmax_30C = nanmean(tempmat(ivec,2:4),2);
grmax_err_30C = nanstd(tempmat(ivec,2:4),0,2);

%% 37C NZCasePlus
tempmat = monod_mat_NZCasePlus_37C;
smat = size(tempmat);


grmax_37C = nanmean(tempmat(ivec,2:4),2);
grmax_err_37C = nanstd(tempmat(ivec,2:4),0,2);

%% Now do Arrhenius analysis
T = [25 30 37]';
clrs = distinguishable_colors(length(ivec));

grmax_tot_NZCasePlus = [grmax_25C grmax_30C grmax_37C];
grmax_err_tot_NZCasePlus = [grmax_err_25C grmax_err_30C grmax_err_37C];

Ea_NZCasePlus = [];
Ea_err_NZCasePlus = [];


cols = cbrewer('div','RdBu', length(concentration));

figure;
for k=1:length(ivec)
    g = grmax_tot_NZCasePlus(k,:)';
    ge = grmax_err_tot_NZCasePlus(k,:)';
 
    col = cols(k,:);
    errorbar(1./(T+273.15), log(g), ge./g, 'o','color', col, 'linewidth', 1);
    hold on;
    
    f = fitlm(1./(T+273.15), log(g),'Weights', 1./(ge./g).^2);
    plot(1./(T+273.15),f.feval(1./(T+273.15)), 'color', col, 'linewidth', 1);
    
    ea = -f.Coefficients.Estimate(2)/298/1.688;
    ea_se = -f.Coefficients.SE(2)/298/1.688;

    disp('Activation Energy (KT):');
    disp(ea);
    
    disp('With error:');
    disp(ea_se);
    
    Ea_NZCasePlus = [Ea_NZCasePlus; ea];
    Ea_err_NZCasePlus = [Ea_err_NZCasePlus; ea_se];
    
    
end

set(gca, 'fontsize', 20);
set(gcf, 'Position', [0 0 400 300]);
xlabel('1/T (1/K)');
ylabel('Log (growth rate)');
xlim([1/312 1/296]);
box off;

%% Plot Ea vs. concentration (in units of kcal/mol)
figure;
errorbar(concentration, Ea_NZCasePlus, Ea_err_NZCasePlus, 'ko-', 'markerfacecolor', 'g');
ylim([0 40]);
set(gca, 'fontsize', 20);
set(gcf, 'Position', [0 0 400 300]);
xlabel('Log(concentration)');
ylabel('Activation energy (kcal/mol)');
box off
% xlim([0 max(concentration)*1.1])

%% Compute average activation energy at saturation, with standard error
sat_index = [1:3]; % Indexing for saturating concentrations

Ea_mean = mean(Ea_NZCasePlus(sat_index))/1.688
Ea_mean_err = sqrt(sum(Ea_err_NZCasePlus(sat_index).^2)/length(sat_index))/1.688
Ea_se = Ea_mean_err/sqrt(length(sat_index))/1.688









